Functions and global parameters
opts_knit$set(root.dir = "..")
# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE, tidy = TRUE)
read_excel_df <- function(...) data.frame(read_excel(...))
trait_space_plot <- function(X, dimensions, group, indices, basecex = 1, title = "", colors = c("#3E4A89FF", "#35B779FF"), point.alpha = 0.7, background.indices = NULL, pch = 1, labels = c("sub-space", "total space"), legend = TRUE) {
total_coors <- as.ppp(as.matrix(X[, dimensions]), c(range(X[, dimensions[1]]), range(X[, dimensions[2]])))
total_space <- raster(density.ppp(total_coors))
xlim <- range(X[, dimensions[1]])
ylim <- range(X[, dimensions[2]])
# ylim <- if (is.null(background.indices))
# range(X[, dimensions[2]]) else
# c(min(X[, dimensions[2]]), c(max(X[, dimensions[2]]) + (max(X[, dimensions[2]]) - min(X[, dimensions[2]])) * 0.11))
plot(x = X[, dimensions[1]], y = X[, dimensions[2]], col = "white", cex = basecex, xlab = dimensions[1], ylab= dimensions[2], xlim = xlim, ylim = ylim, cex.lab = basecex, xaxs="i", yaxs="i")
# add background group
if (!is.null(background.indices)){
# get points
Y_bg <- X[background.indices, ]
if (nrow(Y_bg) > 1 & point.alpha > 0)
points(Y_bg[, dimensions[1]], Y_bg[, dimensions[2]], col = adjustcolor(colors[2], alpha.f = point.alpha), pch = pch)
if (nrow(Y_bg) > 4){
coors <- as.ppp(as.matrix(Y_bg[, dimensions]), c(range(Y_bg[, dimensions[1]]), range(Y_bg[, dimensions[2]])))
raster_dens <- raster(density.ppp(coors))
palpre <- colorRampPalette(c("white", colors[2]))
colsn <- 10
palpre2 <- sapply(1:colsn, function(x) adjustcolor(col = palpre(colsn)[x], alpha = seq(0.1, 0.9, length.out = 20)[x]))
image(raster_dens, add = TRUE, col = palpre2)
}
}
Y <- X[indices, ]
if (nrow(Y) > 1 & point.alpha > 0)
points(Y[, dimensions[1]], Y[, dimensions[2]], col = adjustcolor(colors[1], alpha.f = point.alpha), pch = pch)
if (nrow(Y) > 4){
coors <- as.ppp(as.matrix(Y[, dimensions]), c(range(Y[, dimensions[1]]), range(Y[, dimensions[2]])))
raster_dens <- raster(density.ppp(coors))
palpre <- colorRampPalette(c("white", colors[1]))
colsn <- 10
palpre2 <- sapply(1: colsn, function(x) adjustcolor(col = palpre(colsn)[x], alpha = seq(0.1, 0.9, length.out = 20)[x]))
image(raster_dens, add = TRUE, col = palpre2)
}
usr <- par()$usr
# legend
if (!is.null(background.indices) & legend)
# legend(y = rep(usr[4] + ((usr[4] - usr[3]) * 0.05), 2), x = c(usr[1] + (usr[2] - usr[1]) * 0.25, usr[1] + (usr[2] - usr[1]) * 0.75), col = colors, legend = labels, pch = c(pch, pch), cex = basecex * 1.5, bty = "n", horiz = TRUE, xjust = 0)
legend("topright", col = colors, legend = labels, pch = c(pch, pch), cex = basecex * 1.5, xjust = 0, bg = adjustcolor("white", alpha.f = 0.5))
title(title, cex.main = basecex * 1.2)
}
# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")
cols <- viridis(10, alpha = 0.7)
summary_brm_model <- function(x){
summ <- summary(x)$fixed
fit <- x$fit
betas <- grep("^b_", names(fit@sim$samples[[1]]), value = TRUE)
hdis <- t(sapply(betas, function(y) hdi(fit@sim$samples[[1]][[y]]))
)
summ_table <- data.frame(summ, hdis, iterations = attr(fit, "stan_args")[[1]]$iter, chains = length(attr(fit, "stan_args")))
summ_table <- summ_table[rownames(summ_table) != "Intercept", c("Estimate", "Rhat", "Bulk_ESS", "l.95..CI", "u.95..CI", "iterations", "chains")]
summ_table <- as.data.frame(summ_table)
summ_table$Rhat <- round(summ_table$Rhat, digits = 3)
summ_table$CI_low <- round(unlist(summ_table$l.95..CI), digits = 3)
summ_table$CI_high <- round(unlist(summ_table$u.95..CI), digits = 3)
summ_table$l.95..CI <- summ_table$u.95..CI <- NULL
out <- lapply(betas, function(y) data.frame(variable = y, value = sort(fit@sim$samples[[1]][[y]], decreasing = FALSE)))
posteriors <- do.call(rbind, out)
posteriors <- posteriors[posteriors$variable != "b_Intercept", ]
gg <- ggplot(data = posteriors, aes(y = variable, x = value, fill = stat(quantile))) +
geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = HDInterval::hdi, vline_linetype = 2) +
theme_classic() +
xlim(c(min(summ_table[ , "CI_low"]), max(summ_table[ , "CI_high"]))) +
geom_vline(xintercept = 0, col = "red") +
scale_fill_manual(values = c("transparent", "lightblue", "transparent"), guide = "none") + ggtitle(x$formula)
summ_table$Rhat <- ifelse(summ_table$Rhat > 1.05, cell_spec(summ_table$Rhat, "html", color ="white", background = "red", bold = TRUE, font_size = 12), cell_spec(summ_table$Rhat, "html"))
signif <- summ_table[,"CI_low"] * summ_table[,"CI_high"] > 0
df1 <- kbl(summ_table, row.names = TRUE, escape = FALSE, format = "html", digits = 3)
df1 <- row_spec(kable_input = df1,row = which(summ_table$CI_low * summ_table$CI_high > 0), background = adjustcolor(cols[9], alpha.f = 0.3))
df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12)
print(df1)
print(gg)
}
# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")
# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd", "freq.median",
"freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
"mindom", "maxdom", "dfrange", "modindx", "meanpeakf")]
sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])
sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year, sep = "-")
sels$date <- as.Date(sels$date, format = "%d-%b-%Y")
Exploratory graph
Physiological parameters
breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count", "D14.Breath.Count",
"D21.Breath.Count", "D28.Breath.Count")])
weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.", "D14.Bird.Weight..g.",
"D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])
cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress", "D14.CORT.Stress",
"D21.CORT.Stress", "D28.CORT.Stress")])
cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline", "D14.CORT.Baseline",
"D21.CORT.Baseline", "D28.CORT.Baseline")])
stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round", "Treatment.Room")],
week = breath.count$ind, breath.count = breath.count$values, weight = weight$values,
cort.stress = cort.stress$values, cort.baseline = cort.baseline$values, stress.response = cort.stress$values -
cort.baseline$values)
# head(stress)
stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."), "[[", 1),
levels = c("D3", "D7", "D14", "D21", "D28"))
stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)
stress$treatment <- factor(stress$Treatment, levels = c("Control", "Medium Stress",
"High Stress"))
# remove 5th week
stress <- stress[stress$week != "D28", ]
stress_l <- lapply(stress$Bird.ID, function(x) {
X <- stress[stress$Bird.ID == x, ]
X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
X$weight <- X$weight - X$weight[X$week == "D3"]
X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week == "D3"]
X$stress.response <- X$stress.response - X$stress.response[X$week == "D3"]
return(X)
})
stress <- do.call(rbind, stress_l)
ggplot(data = stress, aes(x = day, y = weight, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Weight (g)", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = breath.count, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Breath count", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = cort.baseline, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Baseline CORT", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = cort.stress, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Stress CORT", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = stress.response, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Stress response", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

# stress_l <- lapply(unique(stress$Bird.ID), function(x){ X <-
# stress[stress$Bird.ID == x, ] X$cort.baseline.change <- X$cort.baseline - X$
# })
standard_error <- function(x) sd(x)/sqrt(length(x))
agg_stress <- aggregate(cort.baseline ~ week + treatment, stress, mean)
agg_stress$se <- aggregate(cort.baseline ~ week + treatment, stress, standard_error)[,
3]
agg_stress$Week <- 1:4
ggplot(data = agg_stress, aes(x = Week, y = cort.baseline, fill = treatment)) + geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = cort.baseline - se, ymax = cort.baseline + se), width = 0.1) +
# geom_line(size = 1.2) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") +
labs(y = "Baseline CORT", x = "Week") + theme_classic(base_size = 30) + theme(legend.position = "none")

FOXP2
foxp2 <- read.csv("./data/raw/stress_IHC_counts.csv")
foxp2$Treatment[grep("HIGH", foxp2$Treatment)] <- "High stress"
foxp2$Treatment[grep("Control", foxp2$Treatment, ignore.case = T)] <- "Control"
foxp2$Treatment[grep("MED", foxp2$Treatment)] <- "Medium stress"
foxp2$Treatment <- factor(foxp2$Treatment, levels = c("Control", "Medium stress",
"High stress"))
ggplot(data = foxp2, aes(x = Treatment, y = FoxP2.Counts.MMSt.Striatum, fill = Treatment)) +
geom_boxplot() + geom_point() + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
labs(y = "% FoxP2 expression", x = "Treatment") + theme_classic(base_size = 30) +
# scale_x_discrete( labels = c(1:4)) +
theme(legend.position = "none")
Behavioral parameters
Call counts per treatment
call_count <- aggregate(selec ~ week + treatment + round, data = sels, length)
call_count$Week <- factor(call_count$week, levels = (unique(call_count$week))[5:1])
call_count$round <- paste("Round", call_count$round)
call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
call_count$stress.response <- aggregate(selec ~ week + treatment + round, data = sels,
mean)[, 4]
ggplot(data = call_count, aes(x = treatment, y = selec, fill = Week)) + geom_bar(stat = "identity") +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + coord_flip() + facet_wrap(~round,
ncol = 2) + labs(y = "Number of calls", x = "Treatment") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.15))

call_count2 <- aggregate(selec ~ week + treatment + round + ID, data = sels, length)
agg_count <- aggregate(selec ~ week + treatment, call_count2, mean)
agg_count$se <- aggregate(selec ~ week + treatment, call_count2, standard_error)[,
3]
agg_count$treatment <- factor(agg_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(data = agg_count, aes(x = week, y = selec, fill = treatment)) + geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = selec - se, ymax = selec + se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Number of calls",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = "none")

Call counts per bird and treatment
# Call counts per treatment by individual (in color)
call_count <- aggregate(selec ~ week + treatment + round + ID, data = sels, length)
call_count$Week <- factor(call_count$week, levels = (1:5))
call_count$round <- paste("Round", call_count$round)
call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(data = call_count, aes(x = week, y = selec, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Call count", x = "Week") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.25))

Acoustic parameters
call_week_params <- aggregate(cbind(duration, meanfreq, sd, freq.median, freq.IQR,
time.IQR, skew, kurt, sp.ent, time.ent, entropy, meandom, mindom, maxdom, dfrange,
modindx, meanpeakf) ~ week + treatment + round, data = sels, mean)
call_week_params_sd <- aggregate(cbind(duration, meanfreq, sd, freq.median, freq.IQR,
time.IQR, skew, kurt, sp.ent, time.ent, entropy, meandom, mindom, maxdom, dfrange,
modindx, meanpeakf) ~ week + treatment + round, data = sels, sd)
names(call_week_params_sd) <- names(call_week_params) <- c("week", "treatment", "round",
"Duration", "Mean_frequency_(kHz)", "SD_frequency_(kHz)", "Median_frequency_(kHz)",
"Interquantile_frequency_range_(kHz)", "Interquantile_time_range_(s)", "Skewness",
"Kurtosis", "Spectral_entropy", "Time_entropy", "Overall entropy", "Mean_dominant_frequency_(kHz)",
"Minimum_dominant_frequency_(kHz)", "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)",
"Modulation_index", "Mean_peak_frequency_(kHz)")
call_week_params_sd <- call_week_params_sd[, c("Duration", "Mean_frequency_(kHz)",
"SD_frequency_(kHz)", "Median_frequency_(kHz)", "Interquantile_frequency_range_(kHz)",
"Interquantile_time_range_(s)", "Skewness", "Kurtosis", "Spectral_entropy", "Time_entropy",
"Overall entropy", "Mean_dominant_frequency_(kHz)", "Minimum_dominant_frequency_(kHz)",
"Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)", "Modulation_index",
"Mean_peak_frequency_(kHz)")]
names(call_week_params_sd) <- paste(names(call_week_params_sd), "sd", sep = ".")
call_week_params <- cbind(call_week_params, call_week_params_sd)
call_week_params$round <- paste("Round", call_week_params$round)
call_week_params$treatment <- factor(call_week_params$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
pd <- position_dodge(0.3)
ggs <- lapply(c("Duration", "Mean_frequency_(kHz)", "SD_frequency_(kHz)", "Median_frequency_(kHz)",
"Interquantile_frequency_range_(kHz)", "Interquantile_time_range_(s)", "Skewness",
"Kurtosis", "Spectral_entropy", "Time_entropy", "Overall entropy", "Mean_dominant_frequency_(kHz)",
"Minimum_dominant_frequency_(kHz)", "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)",
"Modulation_index", "Mean_peak_frequency_(kHz)"), function(i) {
call_week_params$var <- call_week_params[, i]
call_week_params$var.min <- call_week_params$var - call_week_params[, paste(i,
"sd", sep = ".")]
call_week_params$var.max <- call_week_params$var + call_week_params[, paste(i,
"sd", sep = ".")]
ggplot(call_week_params, aes(x = week, y = var, col = treatment)) + geom_point(size = 4,
position = pd) + geom_errorbar(aes(ymin = var.min, ymax = var.max), width = 0.3,
position = pd, size = 1.7) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = i, x = "Week") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.25))
})
ggs
## [[1]]

##
## [[2]]

##
## [[3]]

##
## [[4]]

##
## [[5]]

##
## [[6]]

##
## [[7]]

##
## [[8]]

##
## [[9]]

##
## [[10]]

##
## [[11]]

##
## [[12]]

##
## [[13]]

##
## [[14]]

##
## [[15]]

##
## [[16]]

##
## [[17]]

Acoustic space projection
scale_acous_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median",
"freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
"mindom", "maxdom", "dfrange", "modindx", "meanpeakf")])
urfmod <- randomForest(x = scale_acous_param, importance = TRUE, proximity = TRUE,
ntree = 10000)
saveRDS(urfmod, "./data/processed/unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")
mds_urf_prox <- cmdscale(1 - urfmod$proximity, k = 2)
saveRDS(mds_urf_prox, "./data/processed/MDS_unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")
Random forest
urfmod <- readRDS("./data/processed/unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")
mds_urf_prox <- readRDS("./data/processed/MDS_unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")
sels$MDS1 <- mds_urf_prox[, 1]
sels$MDS2 <- mds_urf_prox[, 2]
print(1 - sum(urfmod$votes[, 1] > urfmod$votes[, 2])/nrow(urfmod$votes))
## [1] 0.005058366
Proportion of real data classified as synthetic data by the unsupervised random forest: 0.005058
ggplot(sels, aes(x = MDS1, y = MDS2, col = as.factor(round))) + geom_point() + labs(color = "Round") +
scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 30) + theme(legend.position = c(0.62,
0.3)) + guides(colour = guide_legend(override.aes = list(size = 10)))

PCA
pca <- prcomp(x = sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
"time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
"maxdom", "dfrange", "modindx", "meanpeakf")], scale. = TRUE)
Y <- as.data.frame(pca$x[, 1:2])
names(Y) <- c("PC1", "PC2")
sels <- data.frame(sels, Y)
ggplot(sels, aes(x = PC1, y = PC2, col = as.factor(round))) + geom_point() + labs(color = "Round") +
scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) + theme(legend.position = c(0.9,
0.8)) + guides(colour = guide_legend(override.aes = list(size = 10)))

t-SNE
scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
"time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
"maxdom", "dfrange", "modindx", "meanpeakf")])
tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE, max_iter = 5000)
saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")
sels <- data.frame(sels, Y)
ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(round))) + geom_point() +
labs(color = "Round") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
theme(legend.position = c(0.955, 0.25)) + guides(colour = guide_legend(override.aes = list(size = 10)))

sels$treatment <- factor(sels$treatment, levels = c("Control", "Medium Stress", "High Stress"))
ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(treatment))) + geom_point() +
labs(color = "Treatment") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
theme(legend.position = c(0.88, 0.2)) + guides(colour = guide_legend(override.aes = list(size = 10)))

Acoustic space by individual (each point) for the 3 dimension reduction methods
- Only individuals with at least 20 calls
- Using rarefaction with n = the minimum sample size across all individuals
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab[tab > 20]), ]
pca_areas <- rarefact_space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID",
parallel = 4, pb = TRUE, type = "density")
rf_mds_areas <- rarefact_space_size(X = Y, dimensions = c("MDS1", "MDS2"), group = "ID",
parallel = 4, pb = TRUE, type = "density")
tsne_areas <- rarefact_space_size(X = Y, dimensions = c("TSNE1", "TSNE2"), group = "ID",
parallel = 4, pb = TRUE, type = "density")
areas <- data.frame(pca_areas[, -ncol(pca_areas)], pca.area = pca_areas$mean.size,
rf.mds.area = rf_mds_areas$mean.size, tsne.area = tsne_areas$mean.size)
saveRDS(areas, "./data/processed/acoustic_space_area_by_individual_for_each_dimension_reduction_method.RDS")
areas <- readRDS("./data/processed/acoustic_space_area_by_individual_for_each_dimension_reduction_method.RDS")
ggpairs(areas, columns = which(names(areas) %in% c("pca.area", "rf.mds.area", "tsne.area")),
columnLabels = c("PCA", "Rand.for. MDS", "t-SNE")) + theme_minimal(base_size = 30)

Individual acoustic space across weeks
- PCA and t-SNE used for dimensionality reduction
- Only weeks with at least 6 calls were included
PCA
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID.week)
# sort(tab)
Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]
# run with rarefaction
areas_by_week_raref <- rarefact_space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
parallel = 1, pb = TRUE, type = "density")
# run without rarefaction
areas_by_week <- space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
parallel = 1, pb = TRUE, type = "density", extent = 2)
areas_by_week$raref.area <- areas_by_week_raref$mean.size
# add metadata
areas_by_week$ID <- sapply(strsplit(areas_by_week$group, "-"), "[[", 1)
areas_by_week$week <- factor(sapply(strsplit(areas_by_week$group, "-"), "[[", 2),
levels = 1:5)
areas_by_week$treatment <- sapply(1:nrow(areas_by_week), function(x) metadat$Treatment[metadat$Bird.ID ==
areas_by_week$ID[x]][1])
areas_by_week$round <- sapply(1:nrow(areas_by_week), function(x) metadat$Round[metadat$Bird.ID ==
areas_by_week$ID[x]][1])
areas_by_week$round <- paste("Round", areas_by_week$round)
call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
areas_by_week$treatment <- factor(areas_by_week$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
saveRDS(areas_by_week, "./data/processed/acoustic_space_area_by_individual_and_week.RDS")
With rarefaction
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

With rarefaction and compared to week 1
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {
X <- areas_by_week[areas_by_week$ID == x, ]
X$raref.area <- X$raref.area - X$raref.area[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction
ggplot(data = areas_by_week[areas_by_week$week != 5, ], aes(x = week, y = size, group = ID,
col = treatment)) + geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction and compared to week 1
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {
X <- areas_by_week[areas_by_week$ID == x, ]
X$size <- X$size - X$size[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
ggplot(data = areas_by_week, aes(x = week, y = size, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Acoustic space area by number of calls per individual and week
ggplot(data = areas_by_week, aes(x = n, y = size, col = treatment)) + geom_point() +
# geom_smooth(size = 1.2, method = 'lm') +
scale_color_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~round, ncol = 2, scales = "free_x") +
labs(y = "Acoustic space area", x = "Number of calls (n)") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.25))

- Acoustic space area by number of calls per individual and week by treatment
- Must have at least 2 calls
ggplot(data = areas_by_week, aes(x = n, y = log(size + 10), col = treatment)) + geom_point() +
geom_smooth(size = 1.2, method = "lm") + scale_color_viridis_d(begin = 0.1, end = 0.9) +
# facet_wrap(~ round, ncol = 2, scales = 'free_x') +
labs(y = "log(acoustic space area)", x = "Number of calls (n)") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.7))

Save acoustic space plots
trait_space_plot(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
"TSNE2"), point.alpha = 0.2, pch = 1)
#
trait_space_plot(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
"TSNE2"), background.indices = which(sels$ID != "395WBHM" & sels$record.group ==
sels$record.group[sels$ID == "395WBHM"][1]), point.alpha = 0.2, pch = 1)
# Y <- sels[sels$ID == '395WBHM',] Y <- sels[sels$ID %in% c('395WBHM',
# '394WBHM', '360YGHM', '385YBHM'), ]
# trait_overlap(Y, dimensions = c('TSNE1', 'TSNE2'), group = 'ID', type =
# 'density', overlap.type = 'assymetric')
# X = Y dimensions = c('TSNE1', 'TSNE2') level <- '395WBHM' basecex <- 1 group
# <- 'ID' indices <- which(X[, group] == level)
trait_space_plot(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
"TSNE2"), background.indices = which(sels$ID != "395WBHM" & sels$record.group ==
sels$record.group[sels$ID == "395WBHM"][1]), point.alpha = 0.2, pch = 1)
# with points
out <- pblapply(unique(sels$ID), function(x) {
tiff(filename = file.path("./images", paste(x, sels$treatment[sels$ID == x][1],
"v2.tiff", sep = "-")), res = 120, width = 800, height = 800)
par(mfrow = c(2, 2))
# W <- sels[sels$ID == x, ]
for (i in 1:4) try(trait_space_plot(X = sels, dimensions = c("TSNE1", "TSNE2"),
indices = which(sels$ID == x & sels$week == i), background.indices = which(sels$ID !=
x & sels$week == i & sels$record.group == sels$record.group[sels$ID ==
x][1]), title = paste(x, "week", i), basecex = 1, labels = c(x, "group"),
legend = FALSE), silent = TRUE)
dev.off()
})
TSNE
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID.week)
# sort(tab)
Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]
# run with rarefaction
areas_by_week_raref <- rarefact_space_size(X = Y, dimensions = c("TSNE1", "TSNE2"),
group = "ID.week", parallel = 1, pb = TRUE, type = "density")
# run without rarefaction
areas_by_week <- space_size(X = Y, dimensions = c("TSNE1", "TSNE2"), group = "ID.week",
parallel = 4, pb = TRUE, type = "density")
areas_by_week$raref.area <- areas_by_week_raref$mean.size
# add metadata
areas_by_week$ID <- sapply(strsplit(areas_by_week$group, "-"), "[[", 1)
areas_by_week$week <- factor(sapply(strsplit(areas_by_week$group, "-"), "[[", 2),
levels = 1:5)
areas_by_week$treatment <- sapply(1:nrow(areas_by_week), function(x) metadat$Treatment[metadat$Bird.ID ==
areas_by_week$ID[x]][1])
areas_by_week$round <- sapply(1:nrow(areas_by_week), function(x) metadat$Round[metadat$Bird.ID ==
areas_by_week$ID[x]][1])
areas_by_week$round <- paste("Round", areas_by_week$round)
call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
areas_by_week$treatment <- factor(areas_by_week$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
saveRDS(areas_by_week, "./data/processed/acoustic_space_area_by_individual_and_week_TSNE.RDS")
With rarefaction
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week_TSNE.RDS")
ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

With rarefaction and compared to week 1
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week_TSNE.RDS")
areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {
X <- areas_by_week[areas_by_week$ID == x, ]
X$raref.area <- X$raref.area - X$raref.area[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction
ggplot(data = areas_by_week[areas_by_week$week != 5, ], aes(x = week, y = size, group = ID,
col = treatment)) + geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction and compared to week 1
areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {
X <- areas_by_week[areas_by_week$ID == x, ]
X$size <- X$size - X$size[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
ggplot(data = areas_by_week, aes(x = week, y = size, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Acoustic space area by number of calls per individual and week
ggplot(data = areas_by_week, aes(x = n, y = size, col = treatment)) + geom_point() +
# geom_smooth(size = 1.2, method = 'lm') +
scale_color_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~round, ncol = 2, scales = "free_x") +
labs(y = "Acoustic space area", x = "Number of calls (n)") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.25))

- Acoustic space area by number of calls per individual and week by treatment
- Must have at least 2 calls
ggplot(data = areas_by_week, aes(x = n, y = log(size + 10), col = treatment)) + geom_point() +
geom_smooth(size = 1.2, method = "lm") + scale_color_viridis_d(begin = 0.1, end = 0.9) +
# facet_wrap(~ round, ncol = 2, scales = 'free_x') +
labs(y = "log(acoustic space area)", x = "Number of calls (n)") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.7))

Acoustic space overlap
- 2 methods: pairwise distance between observations and kernel density overlap
Individual overlap to its own group
PCA
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 2], ]
group_ovlp_l <- lapply(unique(Y$ID.week), function(x) {
# group data
X <- Y[Y$record.group == Y$record.group[Y$ID.week == x][1] & Y$week == Y$week[Y$ID.week ==
x][1], ]
# label all other individuals as group
X$ID2 <- ifelse(X$ID.week == x, "focal", "group")
# run with density
ovlp <- if (min(table(X$ID2)) > 4)
rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"), group = "ID2",
parallel = 1, pb = FALSE, outliers = 0.95, type = "mean.density.overlap") else matrix(rep(NA, 3), nrow = 1)
dsts <- if (min(table(X$ID2)) > 1)
rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"), group = "ID2",
parallel = 1, pb = FALSE, outliers = 0.95, type = "distance") else matrix(rep(NA, 3), nrow = 1)
out <- data.frame(ID = Y$ID[Y$ID.week == x][1], group = Y$record.group[Y$ID.week ==
x][1], week = Y$week[Y$ID.week == x][1], overlap.to.group = ovlp[1, 3], distance.to.group = dsts[1,
3], treatment = Y$treatment[Y$ID.week == x][1], round = Y$round[Y$ID.week ==
x][1])
return(out)
})
group_ovlp <- do.call(rbind, group_ovlp_l)
saveRDS(group_ovlp, "./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")
group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")
group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID), function(x) {
X <- group_ovlp[group_ovlp$ID == x, ]
X$overlap.to.group <- X$overlap.to.group - X$overlap.to.group[X$week == min(as.numeric(as.character(X$week)))]
X$distance.to.group <- X$distance.to.group - X$distance.to.group[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to group acoustic space",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))

ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to group acoustic space",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID.week)
Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]
# run with rarefaction
overlap_by_week <- rarefact_space_similarity(X = Y, dimensions = c("PC1", "PC2"),
group = "ID.week", parallel = 1, pb = TRUE, outliers = 0.95, type = "mean.density.overlap")
saveRDS(overlap_by_week, "./data/processed/acoustic_space_overlap_by_individual_and_week.RDS")
TSNE
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 6], ]
group_ovlp_l <- lapply(unique(Y$ID.week), function(x) {
# group data
X <- Y[Y$record.group == Y$record.group[Y$ID.week == x][1] & Y$week == Y$week[Y$ID.week ==
x][1], ]
# label all other individuals as group
X$ID2 <- ifelse(X$ID.week == x, "focal", "group")
# run with density
ovlp <- if (min(table(X$ID2)) >= 6)
rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"), group = "ID2",
parallel = 1, pb = FALSE, type = "mean.density.overlap")[, 1:3] else matrix(rep(NA, 3), nrow = 1)
dsts <- if (min(table(X$ID2)) >= 2)
rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"), group = "ID2",
parallel = 1, pb = FALSE, type = "distance")[, 1:3] else matrix(rep(NA, 3), nrow = 1)
out <- data.frame(ID = X$ID[X$ID.week == x][1], group = X$record.group[X$ID.week ==
x][1], week = X$week[X$ID.week == x][1], overlap.to.group = ovlp[1, 3], distance.to.group = dsts[1,
3], treatment = X$treatment[X$ID.week == x][1], round = X$round[X$ID.week ==
x][1])
return(out)
})
group_ovlp <- do.call(rbind, group_ovlp_l)
saveRDS(group_ovlp, "./data/processed/acoustic_space_density_overlap_to_group_by_week_TSNE.RDS")
group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week_TSNE.RDS")
# comparing to week 1 group_ovlp <- do.call(rbind,
# lapply(unique(group_ovlp$ID), function(x){ X <- group_ovlp[group_ovlp$ID ==
# x, ] X$overlap.to.group <- X$overlap.to.group / X$overlap.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] X$distance.to.group <-
# X$distance.to.group - X$distance.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] return(X) }))
group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to group acoustic space",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))

Change in overlap with self over time
PCA
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 6], ]
indiv_ovlp_l <- lapply(unique(Y$ID), function(x) {
# group data
X <- Y[Y$ID == x, ]
tab_week <- table(X$week)
X <- X[X$week %in% names(tab_week)[tab_week >= 6], ]
if (any(X$week == 1)) {
# run with density
ovlp <- try(rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"),
group = "week", parallel = 1, pb = FALSE, type = "mean.density.overlap"),
silent = TRUE)
dsts <- try(rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"),
group = "week", parallel = 1, pb = FALSE, type = "distance"), silent = TRUE)
out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = if (is(ovlp,
"try-error"))
NA else ovlp$group.2[ovlp$group.1 == 1], overlap.to.first.week = if (is(ovlp,
"try-error"))
NA else ovlp$mean.overlap[ovlp$group.1 == 1], distance.to.first.week = if (is(dsts,
"try-error"))
NA else dsts$mean.distance[dsts$group.1 == 1], treatment = X$treatment[X$ID ==
x][1], round = X$round[X$ID == x][1])
} else out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = NA,
overlap.to.first.week = NA, distance.to.first.week = NA, treatment = X$treatment[X$ID ==
x][1], round = X$round[X$ID == x][1])
return(out)
})
indiv_ovlp <- do.call(rbind, indiv_ovlp_l)
saveRDS(indiv_ovlp, "./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")
indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")
indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
indiv_ovlp <- indiv_ovlp[!is.na(indiv_ovlp$group), ]
ggplot(data = indiv_ovlp, aes(x = week, y = overlap.to.first.week, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to acoustic space from first week",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))

ggplot(data = indiv_ovlp, aes(x = week, y = distance.to.first.week, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to acoustic space from first week",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))

agg_dist <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, mean)
agg_dist$se <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, standard_error)[,
3]
agg_dist$Week <- 1:4
# ggplot(data = agg_dist, aes(x = Week, y = distance.to.first.week, fill =
# treatment)) + geom_bar(stat = 'identity') +
# geom_errorbar(aes(ymin=distance.to.first.week - se, ymax =
# distance.to.first.week + se), width=.1) + # geom_line(size = 1.2) +
# scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~ treatment, ncol =
# 3, scale = 'fixed') + labs(y = 'Distance to self first week', x = 'Week') +
# theme_classic(base_size = 30) + # scale_x_discrete( labels = c(1:4)) +
# theme(legend.position = 'none')
agg_ovlp <- aggregate(overlap.to.first.week ~ week + treatment, indiv_ovlp, mean)
agg_ovlp$se <- aggregate(overlap.to.first.week ~ week + treatment, indiv_ovlp, standard_error)[,
3]
agg_ovlp$Week <- 1:4
# ggplot(data = agg_ovlp, aes(x = Week, y = overlap.to.first.week, fill =
# treatment)) + geom_bar(stat = 'identity') +
# geom_errorbar(aes(ymin=overlap.to.first.week - se, ymax =
# overlap.to.first.week + se), width=.1) + # geom_line(size = 1.2) +
# scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~ treatment, ncol =
# 3, scale = 'fixed') + labs(y = 'Overlap to \nself first week', x = 'Week') +
# theme_classic(base_size = 30) + # scale_x_discrete( labels = c(1:4)) +
# theme(legend.position = 'none')
TSNE
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 2], ]
indiv_ovlp_l <- lapply(unique(Y$ID), function(x) {
# group data
X <- Y[Y$ID == x, ]
tab_week <- table(X$week)
X <- X[X$week %in% names(tab_week)[tab_week >= 6], ]
if (any(X$week == 1)) {
# run with density
ovlp <- try(rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"),
group = "week", parallel = 1, pb = FALSE, type = "mean.density.overlap"),
silent = TRUE)
dsts <- try(rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"),
group = "week", parallel = 1, pb = FALSE, type = "distance"), silent = TRUE)
out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = if (is(ovlp,
"try-error"))
NA else ovlp$group.2[ovlp$group.1 == 1], overlap.to.first.week = if (is(ovlp,
"try-error"))
NA else ovlp$mean.overlap[ovlp$group.1 == 1], distance.to.first.week = if (is(dsts,
"try-error"))
NA else dsts$mean.distance[dsts$group.1 == 1], treatment = X$treatment[X$ID ==
x][1], round = X$round[X$ID == x][1])
} else out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = NA,
overlap.to.first.week = NA, distance.to.first.week = NA, treatment = X$treatment[X$ID ==
x][1], round = X$round[X$ID == x][1])
return(out)
})
indiv_ovlp <- do.call(rbind, indiv_ovlp_l)
saveRDS(indiv_ovlp, "./data/processed/acoustic_space_density_overlap_to_first_week_by_individual_TSNE.RDS")
indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual_TSNE.RDS")
indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
indiv_ovlp <- indiv_ovlp[!is.na(indiv_ovlp$group), ]
ggplot(data = indiv_ovlp, aes(x = week, y = overlap.to.first.week, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to acoustic space from first week",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))

ggplot(data = indiv_ovlp, aes(x = week, y = distance.to.first.week, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to acoustic space from first week",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))

agg_dist <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, mean)
agg_dist$se <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, standard_error)[,
3]
# agg_stress$Week <- 1:4
# ggplot(data = agg_dist, aes(x = week, y = distance.to.first.week, fill =
# treatment)) + geom_bar(stat = 'identity') +
# geom_errorbar(aes(ymin=distance.to.first.week - se, ymax =
# distance.to.first.week + se), width=.1) + # geom_line(size = 1.2) +
# scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~ treatment, ncol =
# 3, scale = 'fixed') + labs(y = 'Distance to first week', x = 'Week') +
# theme_classic(base_size = 30) + # scale_x_discrete( labels = c(1:4)) +
# theme(legend.position = 'none')
Proposed Models for Stress and Learning Pilot Data based on 9/28/21 discussion
Model set 1: effects of treatment
Predicted variables: - Physiology: Baseline CORT, stress response CORT (stress), difference baseline-stress response CORT (stress response), weight, breath rate - Learning: number of calls, total acoustic area, distance moved from week 1, overlap to self week 1, overlap self to group
General form: - Predicted ~ treatment + week + ind(rand) - Predicted ~ treatment + week + treat * week + ind(rand) - Predicted ~ round + week [to check for effect of round]
Note: Use week 1 as baseline, comparing weeks 2, 3, 4
agg_dat <- aggregate(selec ~ ID + week, data = sels, length)
# compare to week 1 agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) {
# baseline <- agg_dat$selec[agg_dat$week == 1 & agg_dat$ID == agg_dat$ID[x]] if
# (length(baseline) > 0) change <- agg_dat$selec[x] - baseline else change <-
# agg_dat$selec[x] return(change) } )
# without comparing to week 1
agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) agg_dat$selec[x])
agg_dat$selec <- NULL
agg_dat$baseline.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.baseline[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.baseline[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.response <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$stress.response[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$stress.response[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.stress[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.stress[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$weight <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$weight[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$weight[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$breath.rate <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$breath.count[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$breath.count[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$acoustic.area <- sapply(1:nrow(agg_dat), function(x) {
area <- areas_by_week$raref.area[areas_by_week$ID == agg_dat$ID[x] & areas_by_week$week ==
agg_dat$week[x]]
if (length(area) < 1)
area <- NA
return(area)
})
agg_dat$acoustic.distance <- sapply(1:nrow(agg_dat), function(x) {
distance <- indiv_ovlp$distance.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
indiv_ovlp$week == agg_dat$week[x]]
if (length(distance) < 1)
distance <- NA
return(distance)
})
agg_dat$acoustic.overlap <- sapply(1:nrow(agg_dat), function(x) {
overlap <- indiv_ovlp$overlap.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
indiv_ovlp$week == agg_dat$week[x]]
if (length(overlap) < 1)
overlap <- NA
return(overlap)
})
agg_dat$acoustic.overlap.to.group <- sapply(1:nrow(agg_dat), function(x) {
overlap <- group_ovlp$overlap.to.group[group_ovlp$ID == agg_dat$ID[x] & group_ovlp$week ==
agg_dat$week[x]]
if (length(overlap) < 1)
overlap <- NA
return(overlap)
})
agg_dat$treatment <- sapply(1:nrow(agg_dat), function(x) unique(sels$treatment[sels$ID ==
agg_dat$ID[x]]))
agg_dat$round <- sapply(1:nrow(agg_dat), function(x) unique(sels$round[sels$ID ==
agg_dat$ID[x]]))
# remove week 1
agg_dat <- agg_dat[agg_dat$week != 1, ]
responses <- setdiff(names(agg_dat), c("ID", "week", "treatment", "round"))
predictors <- c("~ treatment + week + (1|ID)", "~ treatment * week + (1|ID)", "~ round + week + (1|ID)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
vars_to_scale <- names(agg_dat)[sapply(agg_dat, is.numeric)]
vars_to_scale <- vars_to_scale[!vars_to_scale %in% c("week", "round")]
for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[, vars_to_scale])
treatment_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
]
sub_dat
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 2), silent = TRUE)
return(mod)
})
names(treatment_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(treatment_models, "./data/processed/treatment_models.RDS")
agg_dat$week <- as.factor(agg_dat$week)
treatment_models_week_factor <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
]
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 2), silent = TRUE)
return(mod)
})
names(treatment_models_week_factor) <- paste(formulas$responses, formulas$predictors)
saveRDS(treatment_models_week_factor, "./data/processed/treatment_models_week_factor.RDS")
Takeaways
- No significant nor marginal effects
Model set 2: individual stress and effects on learning
General form: - Learning variables ~ baseline CORT [Week 2 immediate stress and learning, week 3 chronic stress and learn]
agg_dat <- agg_dat[agg_dat$week != "4", ]
agg_dat$week <- as.factor(agg_dat$week)
responses <- c("call.count", "acoustic.area", "acoustic.distance", "acoustic.overlap",
"acoustic.overlap.to.group")
predictors <- c("~ baseline.CORT + week + (1|ID)", "~ stress.response + week + (1|ID)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
treatment_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
]
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 1), silent = TRUE)
return(mod)
})
names(treatment_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(treatment_models, "./data/processed/stress_on_learning_models.RDS")
treatment_models <- readRDS("./data/processed/stress_on_learning_models.RDS")
for (x in treatment_models) {
summary_brm_model(x)
}
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
-0.192
|
1.002
|
997.111
|
5000
|
1
|
-0.546
|
0.177
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
-0.154
|
1
|
537.359
|
5000
|
1
|
-0.462
|
0.159
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
0.007
|
1.015
|
256.643
|
5000
|
1
|
-0.248
|
0.274
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
0.101
|
1
|
616.351
|
5000
|
1
|
-0.247
|
0.44
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
0.299
|
1.001
|
1394.371
|
5000
|
1
|
-0.112
|
0.69
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
0.013
|
1
|
556.089
|
5000
|
1
|
-0.25
|
0.287
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
-0.058
|
1.001
|
677.164
|
5000
|
1
|
-0.391
|
0.286
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
-0.138
|
1.002
|
869.593
|
5000
|
1
|
-0.49
|
0.209
|

Model set 3: individual stress and effects on learning on week 2 and 3 independently
Week 2 immediate stress and learning, week 3 chronic stress and learn
General form: - Learning variables ~ baseline CORT - Learning variables ~ difference baseline-stress CORT
Week 2
agg_dat_w2 <- agg_dat[agg_dat$week == "2", ]
responses <- c("call.count", "acoustic.area", "acoustic.distance", "acoustic.overlap",
"acoustic.overlap.to.group")
predictors <- c("~ baseline.CORT + (1|ID)", "~ stress.response + (1|ID)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
treatment_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat_w2[!is.na(agg_dat_w2[names(agg_dat_w2) == formulas$responses[x]]),
]
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 1), silent = TRUE)
return(mod)
})
names(treatment_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(treatment_models, "./data/processed/stress_on_learning_models_week_2.RDS")
treatment_models <- readRDS("./data/processed/stress_on_learning_models_week_2.RDS")
for (x in treatment_models) {
summary_brm_model(x)
}
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
-0.177
|
1
|
1201.316
|
5000
|
1
|
-0.501
|
0.164
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
0.037
|
1.002
|
701.318
|
5000
|
1
|
-0.328
|
0.408
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
-0.138
|
1.012
|
135.412
|
5000
|
1
|
-0.426
|
0.133
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
-0.03
|
1.045
|
78.7
|
5000
|
1
|
-0.329
|
0.298
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
0.11
|
1.004
|
1054.756
|
5000
|
1
|
-0.257
|
0.458
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
0.261
|
1.001
|
831.794
|
5000
|
1
|
-0.072
|
0.586
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
-0.178
|
1
|
227.403
|
5000
|
1
|
-0.545
|
0.202
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
0.013
|
1
|
1041.844
|
5000
|
1
|
-0.232
|
0.25
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
-0.045
|
1.002
|
546.832
|
5000
|
1
|
-0.381
|
0.267
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
-0.138
|
1.004
|
1319.857
|
5000
|
1
|
-0.529
|
0.261
|

Week 3
agg_dat_w3 <- agg_dat[agg_dat$week == "3", ]
responses <- c("call.count", "acoustic.area", "acoustic.distance", "acoustic.overlap",
"acoustic.overlap.to.group")
predictors <- c("~ baseline.CORT + (1|ID)", "~ stress.response + (1|ID)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
treatment_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat_w3[!is.na(agg_dat_w3[names(agg_dat_w3) == formulas$responses[x]]),
]
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 1), silent = TRUE)
return(mod)
})
names(treatment_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(treatment_models, "./data/processed/stress_on_learning_models_week_3.RDS")
treatment_models <- readRDS("./data/processed/stress_on_learning_models_week_3.RDS")
for (x in treatment_models) {
summary_brm_model(x)
}
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
-0.192
|
1.002
|
997.111
|
5000
|
1
|
-0.546
|
0.177
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
-0.154
|
1
|
537.359
|
5000
|
1
|
-0.462
|
0.159
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
0.007
|
1.015
|
256.643
|
5000
|
1
|
-0.248
|
0.274
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
baseline.CORT
|
0.101
|
1
|
616.351
|
5000
|
1
|
-0.247
|
0.44
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
0.299
|
1.001
|
1394.371
|
5000
|
1
|
-0.112
|
0.69
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
0.013
|
1
|
556.089
|
5000
|
1
|
-0.25
|
0.287
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
-0.058
|
1.001
|
677.164
|
5000
|
1
|
-0.391
|
0.286
|
|
|
Estimate
|
Rhat
|
Bulk_ESS
|
iterations
|
chains
|
CI_low
|
CI_high
|
|
stress.response
|
-0.138
|
1.002
|
869.593
|
5000
|
1
|
-0.49
|
0.209
|

Takeaways
- No significant nor marginal effects
Session information
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=pt_BR.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] PhenotypeSpace_0.1.0 warbleR_1.1.26 NatureSounds_1.0.4
## [4] seewave_2.1.8 tuneR_1.3.3.1 ggridges_0.5.3
## [7] posterior_1.2.0 ggrepel_0.9.1 cowplot_1.1.1
## [10] tidybayes_3.0.2 modelr_0.1.8 tidyr_1.2.0
## [13] forcats_0.5.1 purrr_0.3.4 dplyr_1.0.8
## [16] lme4_1.1-28 Matrix_1.3-4 brms_2.16.3
## [19] Rcpp_1.0.8 vegan_2.5-7 lattice_0.20-44
## [22] permute_0.9-7 raster_3.5-15 spatstat_2.3-3
## [25] spatstat.linnet_2.3-2 spatstat.core_2.4-0 rpart_4.1-15
## [28] nlme_3.1-152 spatstat.random_2.1-0 spatstat.geom_2.3-2
## [31] spatstat.data_2.1-2 GGally_2.1.2 adehabitatHR_0.4.19
## [34] adehabitatLT_0.3.25 CircStats_0.2-6 boot_1.3-28
## [37] adehabitatMA_0.3.14 ade4_1.7-18 deldir_1.0-6
## [40] sp_1.4-6 MASS_7.3-54 Rtsne_0.15
## [43] randomForest_4.7-1 formatR_1.11 knitr_1.37
## [46] kableExtra_1.3.4 ggplot2_3.3.5 viridis_0.6.2
## [49] viridisLite_0.4.0 pbapply_1.5-0 readxl_1.3.1
## [52] lubridate_1.8.0 remotes_2.4.2
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 tidyselect_1.1.2 fftw_1.0-6.1
## [4] htmlwidgets_1.5.4 grid_4.1.0 munsell_0.5.0
## [7] codetools_0.2-18 DT_0.20 miniUI_0.1.1.1
## [10] withr_2.4.3 Brobdingnag_1.2-7 colorspace_2.0-3
## [13] highr_0.9 rstudioapi_0.13 stats4_4.1.0
## [16] dtw_1.22-3 tensor_1.5 bayesplot_1.8.1
## [19] labeling_0.4.2 rstan_2.21.3 polyclip_1.10-0
## [22] farver_2.1.0 bridgesampling_1.1-2 coda_0.19-4
## [25] vctrs_0.3.8 generics_0.1.2 xfun_0.29
## [28] R6_2.5.1 markdown_1.1 HDInterval_0.2.2
## [31] gamm4_0.2-6 projpred_2.0.2 bitops_1.0-7
## [34] spatstat.utils_2.3-0 reshape_0.8.8 assertthat_0.2.1
## [37] promises_1.2.0.1 scales_1.1.1 gtable_0.3.0
## [40] processx_3.5.2 goftest_1.2-3 rlang_1.0.1
## [43] systemfonts_1.0.4 splines_4.1.0 broom_0.7.12
## [46] checkmate_2.0.0 inline_0.3.19 yaml_2.3.5
## [49] reshape2_1.4.4 abind_1.4-5 threejs_0.3.3
## [52] crosstalk_1.2.0 backports_1.4.1 httpuv_1.6.5
## [55] rsconnect_0.8.25 tensorA_0.36.2 tools_4.1.0
## [58] ellipsis_0.3.2 jquerylib_0.1.4 RColorBrewer_1.1-2
## [61] proxy_0.4-26 plyr_1.8.6 base64enc_0.1-3
## [64] RCurl_1.98-1.6 ps_1.6.0 prettyunits_1.1.1
## [67] zoo_1.8-9 cluster_2.1.2 magrittr_2.0.2
## [70] ggdist_3.1.0 colourpicker_1.1.1 mvtnorm_1.1-3
## [73] matrixStats_0.61.0 shinyjs_2.1.0 mime_0.12
## [76] evaluate_0.15 arrayhelpers_1.1-0 xtable_1.8-4
## [79] shinystan_2.5.0 gridExtra_2.3 rstantools_2.1.1
## [82] compiler_4.1.0 tibble_3.1.6 crayon_1.5.0
## [85] minqa_1.2.4 StanHeaders_2.21.0-7 htmltools_0.5.2
## [88] mgcv_1.8-36 later_1.3.0 RcppParallel_5.1.5
## [91] DBI_1.1.1 cli_3.2.0 parallel_4.1.0
## [94] igraph_1.2.11 pkgconfig_2.0.3 signal_0.7-7
## [97] terra_1.5-21 spatstat.sparse_2.1-0 xml2_1.3.3
## [100] svUnit_1.0.6 dygraphs_1.1.1.6 svglite_2.1.0
## [103] bslib_0.3.1 webshot_0.5.2 rvest_1.0.2
## [106] stringr_1.4.0 distributional_0.3.0 callr_3.7.0
## [109] digest_0.6.29 rmarkdown_2.11 cellranger_1.1.0
## [112] shiny_1.7.1 gtools_3.9.2 rjson_0.2.21
## [115] nloptr_2.0.0 lifecycle_1.0.1 jsonlite_1.8.0
## [118] fansi_1.0.2 pillar_1.7.0 loo_2.4.1
## [121] fastmap_1.1.0 httr_1.4.2 pkgbuild_1.3.1
## [124] glue_1.6.1 xts_0.12.1 shinythemes_1.2.0
## [127] stringi_1.7.6 sass_0.4.0